---
title: "Ant Swarm Models"
author: "Laura Gomez-Murillo"

---
# This document contains all models done for the paper "Traits that structure birds following army ant swarms and the roles of both facilitation and competition"

#All code ran in R version 4.1.3 

#Preparing R and data:
#1) setting a working directory and uploading libraries and data

library(ggplot2)
library(dplyr)
library(tidyr)
library(tidyverse)
library(vegan)
library(glmmTMB)
library(MuMIn)
library(DHARMa)
library(cowplot)
library(lme4)
library(ggplot2)
library(knitr)
library(DHARMa)

traits <-read.csv("GomezMurillo_2026_Antswarmmech.csv")

# Scaling ant swarm traits and formatting other variables and droping rows that have NA values

s.traits <- traits %>% 
  mutate(s.biomass= scale(biomass)) %>% 
  mutate(s.velocity= scale(velocity)) %>% 
  mutate(ant= as.factor(ant)) %>%
  mutate(year= as.factor(year)) %>%
  mutate(s.total_area= scale(total_area)) %>%
  mutate(s.DominantsAbundance= scale(DominantsAbundance)) %>% 
  mutate(site=as.factor(site)) %>% 
  mutate(swarm_id=as.factor(swarm_id)) %>% 
  mutate(month=as.factor(month)) %>% #made this month rather than date   mutate(ant=as.factor(ant)) %>% 
  mutate(DominantBirds=as.factor(DominantBirds))%>% 
  mutate(pielou.ev=replace(pielou.ev, pielou.ev==0,0.0000000001)) %>% 
  mutate(pielou.ev=replace(pielou.ev, pielou.ev==1,0.9999999999)) %>% 
  mutate(p_n.species=replace(p_n.species, p_n.species==0,0.0000000001)) %>% 
  mutate(p_n.species=replace(p_n.species, p_n.species==1,0.9999999999)) %>%
  mutate(p_shannon=replace(p_shannon, p_shannon==0,0.0000000001)) %>% 
  mutate(p_shannon=replace(p_shannon, p_shannon==1,0.9999999999)) %>% 
  mutate(shannon.di=replace(shannon.di, shannon.di==0,0.0000000001)) %>% 
  mutate(shannon.di=replace(shannon.di, shannon.di==1,0.9999999999)) %>% 
  drop_na()

# Since we are modeling each ant species apart, we have to split data set by ant species
E.bur <- s.traits %>% 
  filter(ant=='Eciton burchellii')
L.pra <- s.traits %>% 
  filter(ant=='Labidus praedator') 

###============ analysis for Eciton burchellii======================#####

### =========== Evenness (Pielou) analysis

#test for different random effect structures
gmPielouEciton=glmmTMB(pielou.ev ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2) + (1|swarm_id), data = E.bur, family=beta_family)
summary(gmPielouEciton)

gmbPielouEciton=glmmTMB(pielou.ev ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2) + (1|site)+(1|swarm_id),data = E.bur, family=beta_family)
summary(gmbPielouEciton)

gmcPielouEciton=glmmTMB(pielou.ev ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2) + (1|site), data = E.bur, family=beta_family)
summary(gmcPielouEciton)

gmdPielouEciton=glmmTMB(pielou.ev ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2) + (1|site)+(1|month), data = E.bur, family=beta_family)
summary(gmdPielouEciton)

gmePielouEciton=glmmTMB(pielou.ev ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2) + (1|swarm_id)+(1|month),data = E.bur, family=beta_family)
summary(gmePielouEciton)

gmfPielouEciton=glmmTMB(pielou.ev ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2) + (1|month),data = E.bur, family=beta_family)
summary(gmfPielouEciton)

#random effect of swarm_id was the best one

#testing residuals of model
e <- qqnorm(resid(gmPielouEciton),pch=16); qqline(resid(gmPielouEciton))

#run model averaging on the model

options(na.action= "na.fail")
mslist <- dredge(gmPielouEciton, subset = dc(s.DominantsAbundance, I(s.DominantsAbundance^2)))
mslist
topmod<-get.models(mslist, subset=cumsum(weight) <= .95)
topmod # here you can see all of the models that make up 95% of the model weight
avgm <- model.avg(topmod)

summary(avgm)
confint(avgm, full = T) 

# predicting for the dominance with the model average

quantile(E.bur$s.DominantsAbundance,c(0.1,0.5,0.90,0.95)) # we graphed to the 95th percentile of the data

s.DominantsAbundance=seq(from=min(E.bur$s.DominantsAbundance),to=quantile(E.bur$s.DominantsAbundance,c(0.95)) ,by=0.1)
s.biomass= mean(E.bur$s.biomass)
s.velocity= mean(E.bur$s.velocity)
s.total_area= mean(E.bur$s.total_area)

newdata <- data.frame(expand.grid(swarm_id=unique(E.bur$swarm_id), s.DominantsAbundance=s.DominantsAbundance,s.biomass=s.biomass, s.velocity=s.velocity, s.total_area=s.total_area))

newdata$fitlink <- predict(avgm, newdata,re.form=NA,se.fit = TRUE, type="link")$fit
newdata$selink <- predict(avgm, newdata, re.form=NA,se.fit = TRUE, type="link")$se.fit
newdata$realval <- plogis(newdata$fitlink)
newdata$upper.ci <- plogis(newdata$fitlink + 1.96*newdata$selink)
newdata$lower.ci <- plogis(newdata$fitlink - 1.96*newdata$selink)
newdata$origlength<-newdata$s.DominantsAbundance * attr(E.bur$s.DominantsAbundance, 'scaled:scale') + attr(E.bur$s.DominantsAbundance, 'scaled:center')
newdata 


#graph the results

pielou_plotting_data <- subset(E.bur, E.bur$DominantsAbundance < quantile(E.bur$DominantsAbundance, 0.95))

pielou <-ggplot(data= newdata, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "Number of dominant individuals", y = "Evenness")+
  scale_x_continuous(expand = expansion(mult = c(0, 0))) + scale_y_continuous(expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#88419d") +
  geom_line() +
  theme_classic()+ 
  theme(text = element_text(size = 18))
pielou
ggsave("PielouplotAug2025.png", plot = pielou, width = 8, height = 6, dpi = 150)


### =========== Species richness model

#test for different random effect structures
E.bur$log_Total_SpeciesSite <- log(E.bur$Total_SpeciesSite)

poisson_Spp=glmmTMB(n_species ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2)  + offset(log_Total_SpeciesSite) + (1|swarm_id),data = E.bur, family = poisson)
summary(poisson_Spp)

poisson_Sppb=glmmTMB(n_species ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2)  + offset(log_Total_SpeciesSite) + (1|site)+(1|swarm_id),data = E.bur, family = poisson)
summary(poisson_Sppb)

poisson_Sppc=glmmTMB(n_species ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2)  + offset(log_Total_SpeciesSite) + (1|site),data = E.bur, family = poisson)
summary(poisson_Sppc)

poisson_Sppd=glmmTMB(n_species ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2)  + offset(log_Total_SpeciesSite) + (1|site)+(1|month),data = E.bur, family = poisson)
summary(poisson_Sppd)

poisson_Sppe=glmmTMB(n_species ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2)  + offset(log_Total_SpeciesSite) + (1|swarm_id)+(1|month),data = E.bur, family = poisson)
summary(poisson_Sppe)

poisson_Sppf=glmmTMB(n_species ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2)  + offset(log_Total_SpeciesSite) + (1|month),data = E.bur, family = poisson)
summary(poisson_Sppf)

#swarm ID is the best one
poisson_Spp=glmmTMB(n_species ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2)  + offset(log_Total_SpeciesSite) + (1|swarm_id),data = E.bur, family = poisson)
summary(poisson_Spp)

confint(poisson_Spp)

#testing residuals of model
s <- qqnorm(resid(poisson_Spp),pch=16); qqline(resid(poisson_Spp))
testDispersion(poisson_Spp)
simulationOutput <- simulateResiduals(fittedModel = poisson_Spp, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

E1 <- resid(poisson_Spp, type = "pearson")
N  <- nrow(E.bur)
p  <- length(fixef(poisson_Spp)) + 1
Overdispersion <- sum(E1^2) / (N - p)
Overdispersion		


#having issues with dredge with the offset. So instead, I am creating my own model list.
poisson_Spp=glmmTMB(n_species ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + I(s.DominantsAbundance^2)  + offset(log_Total_SpeciesSite) + (1|swarm_id),data = E.bur, family = poisson)

options(na.action = "na.fail")

# Data and random structure
data_used <- E.bur
random_effect <- "(1 | swarm_id)"
offset_term <- "offset(log_Total_SpeciesSite)"

# Fixed effects components
main_terms <- c("s.velocity", "s.biomass", "s.total_area", "s.DominantsAbundance")
interaction_terms <- c("s.biomass:s.total_area")
quad_term <- "I(s.DominantsAbundance^2)"

# Create all combinations of main effects
main_combos <- unlist(
  lapply(1:length(main_terms), function(i) combn(main_terms, i, simplify = FALSE)),
  recursive = FALSE
)

# Add optional interaction and quadratic terms
model_formulas <- list()

for (combo in main_combos) {
  # Add interaction only if both s.biomass and s.total_area are in combo
  allow_interaction <- all(c("s.biomass", "s.total_area") %in% combo)

  interaction_add <- if (allow_interaction) c("", interaction_terms) else ""

  # Add quadratic only if s.DominantsAbundance is in combo
  allow_quad <- "s.DominantsAbundance" %in% combo
  quad_add <- if (allow_quad) c("", quad_term) else ""

  # Combine into full model formulas
  for (int_term in interaction_add) {
    for (quad in quad_add) {
      terms <- c(combo, if (int_term != "") int_term, if (quad != "") quad)
      rhs <- paste(c(terms, offset_term, random_effect), collapse = " + ")
      formula_str <- paste("n_species ~", rhs)
      model_formulas[[length(model_formulas) + 1]] <- as.formula(formula_str)
    }
  }
}

model_list <- list()
for (i in seq_along(model_formulas)) {
  f <- model_formulas[[i]]
  fit <- try(glmmTMB(f, data = data_used, family = poisson), silent = TRUE)
  model_list[[i]] <- fit
}
valid_models <- model_list[sapply(model_list, function(m) inherits(m, "glmmTMB"))]
model_selection_table <- model.sel(valid_models)

top_models_95 <- get.models(model_selection_table, subset = cumsum(model_selection_table$weight) <= 0.95)
avg_model <- model.avg(top_models_95)
summary(avg_model)

top_models_95
model_selection_table

#Because no top models emerged (AICc weight > 90%), we performed model averaging 
summary(avg_model)
confint(avg_model, full = T)

#Predicting for the dominance with the model average

scale_center <- attr(E.bur$s.DominantsAbundance, "scaled:center")
scale_scale <- attr(E.bur$s.DominantsAbundance, "scaled:scale")

# Create the scaled variable 
s.DominantsAbundance <- seq(
  from = min(E.bur$s.DominantsAbundance),
  to = quantile(E.bur$s.DominantsAbundance, 0.95),
  by = 0.1
)

s.DominantsAbundance=seq(from=min(E.bur$s.DominantsAbundance),to=quantile(E.bur$s.DominantsAbundance,0.95),by=0.1)
s.biomass= mean(E.bur$s.biomass)
s.velocity= mean(E.bur$s.velocity)
s.total_area= mean(E.bur$s.total_area)
log_Total_SpeciesSite <- mean(E.bur$log_Total_SpeciesSite)

newdata <- data.frame(expand.grid(swarm_id=unique(E.bur$swarm_id), s.DominantsAbundance=s.DominantsAbundance,s.biomass=s.biomass, s.velocity=s.velocity, s.total_area=s.total_area,log_Total_SpeciesSite=log_Total_SpeciesSite))

newdata$fitlink <- predict(avg_model, newdata,re.form=NA,se.fit = TRUE, type="link")$fit
newdata$selink <- predict(avg_model, newdata, re.form=NA,se.fit = TRUE, type="link")$se.fit
newdata$realval <- exp(newdata$fitlink)
newdata$upper.ci <- exp(newdata$fitlink + 1.96*newdata$selink)
newdata$lower.ci <- exp(newdata$fitlink - 1.96*newdata$selink)
newdata$origlength<-newdata$s.DominantsAbundance * attr(E.bur$s.DominantsAbundance, 'scaled:scale') + attr(E.bur$s.DominantsAbundance, 'scaled:center')
newdata 

#plot
spp <- ggplot(data= newdata, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "", y = "Species Richness")+
  scale_x_continuous(expand = expansion(mult = c(0, 0))) + scale_y_continuous(expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#88419d") +
  geom_line() +
  theme_classic()+
  theme(axis.text.x=element_blank())+ # label left in blank to stack all graphs later
  theme(text = element_text(size = 18))   
spp
ggsave("SpeciesRichnessplotSept2025.png", plot = spp, width = 8, height = 6, dpi = 150)

###============== Shannon diversity index ===============####3

#test for different random effect structures
E.bur$log_max_ShannonSite <- log(E.bur$max_ShannonSite)
Gaussian_Shannon <- glmmTMB(
  shannon.di ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + 
    I(s.DominantsAbundance^2) + offset(log_max_ShannonSite) + (1 | swarm_id),
  data = E.bur, family = gaussian)
summary(Gaussian_Shannon)

Gaussian_Shannonb <- glmmTMB(
  shannon.di ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + 
    I(s.DominantsAbundance^2) + offset(log_max_ShannonSite)+ (1|site) + (1 | swarm_id),
  data = E.bur, family = gaussian)
summary(Gaussian_Shannonb)

Gaussian_Shannonc <- glmmTMB(
  shannon.di ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + 
    I(s.DominantsAbundance^2) + offset(log_max_ShannonSite)+ (1|site),
  data = E.bur, family = gaussian)
summary(Gaussian_Shannonc)

Gaussian_Shannond <- glmmTMB(
  shannon.di ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + 
    I(s.DominantsAbundance^2) + offset(log_max_ShannonSite)+ (1|site)+(1|month),
  data = E.bur, family = gaussian)
summary(Gaussian_Shannond)
  
Gaussian_Shannone <- glmmTMB(
shannon.di ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + 
I(s.DominantsAbundance^2) + offset(log_max_ShannonSite)+ (1|swarm_id)+(1|month),
data = E.bur, family = gaussian)
summary(Gaussian_Shannone)

Gaussian_Shannonf <- glmmTMB(
shannon.di ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + 
I(s.DominantsAbundance^2) + offset(log_max_ShannonSite)+ (1|month),
data = E.bur, family = gaussian)
summary(Gaussian_Shannonf)

#swarm ID is the best one
Gaussian_Shannon <- glmmTMB(
  shannon.di ~ s.velocity + s.biomass * s.total_area + s.DominantsAbundance + 
    I(s.DominantsAbundance^2) + offset(log_max_ShannonSite) + (1 | swarm_id),
  data = E.bur, family = gaussian)

summary(Gaussian_Shannon)
confint(Gaussian_Shannon)

#testing residuals of model
sdi<-qqnorm(resid(Gaussian_Shannon),pch=16); qqline(resid(Gaussian_Shannon))
testDispersion(Gaussian_Shannon)
simulationOutput <- simulateResiduals(fittedModel = Gaussian_Shannon, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

# model selection

# Always set this for MuMIn compatibility
options(na.action = "na.fail")

# Data and random structure
data_used <- E.bur
random_effect <- "(1 | swarm_id)"
offset_term <- "offset(log_max_ShannonSite)"

# Fixed effects components
main_terms <- c("s.velocity", "s.biomass", "s.total_area", "s.DominantsAbundance")
interaction_terms <- c("s.biomass:s.total_area")
quad_term <- "I(s.DominantsAbundance^2)"

# Create all combinations of main effects
main_combos <- unlist(
  lapply(1:length(main_terms), function(i) combn(main_terms, i, simplify = FALSE)),
  recursive = FALSE
)

# Add optional interaction and quadratic terms
model_formulas <- list()

for (combo in main_combos) {
  # Add interaction only if both s.biomass and s.total_area are in combo
  allow_interaction <- all(c("s.biomass", "s.total_area") %in% combo)

  interaction_add <- if (allow_interaction) c("", interaction_terms) else ""

  # Add quadratic only if s.DominantsAbundance is in combo
  allow_quad <- "s.DominantsAbundance" %in% combo
  quad_add <- if (allow_quad) c("", quad_term) else ""

  # Combine into full model formulas
  for (int_term in interaction_add) {
    for (quad in quad_add) {
      terms <- c(combo, if (int_term != "") int_term, if (quad != "") quad)
      rhs <- paste(c(terms, offset_term, random_effect), collapse = " + ")
      formula_str <- paste("n_species ~", rhs)
      model_formulas[[length(model_formulas) + 1]] <- as.formula(formula_str)
    }
  }
}

model_list <- list()
for (i in seq_along(model_formulas)) {
  f <- model_formulas[[i]]
  fit <- try(glmmTMB(f, data = data_used, family = poisson), silent = TRUE)
  model_list[[i]] <- fit
}
valid_models <- model_list[sapply(model_list, function(m) inherits(m, "glmmTMB"))]
model_selection_table <- model.sel(valid_models)

top_models_95 <- get.models(model_selection_table, subset = cumsum(model_selection_table$weight) <= 0.95)
avg_model <- model.avg(top_models_95)
summary(avg_model)

top_models_95
model_selection_table

# Because no top models emerged (AICc weight > 90%), we performed model averaging 
summary(avg_model)
confint(avg_model, full = T) 

#Predicting for the dominance with the model average

scale_center <- attr(E.bur$s.DominantsAbundance, "scaled:center")
scale_scale <- attr(E.bur$s.DominantsAbundance, "scaled:scale")

# Create the scaled variable correctly
s.DominantsAbundance <- seq(
  from = min(E.bur$s.DominantsAbundance),
  to = quantile(E.bur$s.DominantsAbundance, 0.95),
  by = 0.1
)

s.DominantsAbundance=seq(from=min(E.bur$s.DominantsAbundance),to=quantile(E.bur$s.DominantsAbundance,0.95),by=0.1)
s.biomass= mean(E.bur$s.biomass)
s.velocity= mean(E.bur$s.velocity)
s.total_area= mean(E.bur$s.total_area)
log_max_ShannonSite <- mean(E.bur$log_max_ShannonSite)

newdata <- data.frame(expand.grid(swarm_id=unique(E.bur$swarm_id), s.DominantsAbundance=s.DominantsAbundance,s.biomass=s.biomass, s.velocity=s.velocity, s.total_area=s.total_area,log_max_ShannonSite=log_max_ShannonSite))

newdata$fitlink <- predict(avg_model, newdata,re.form=NA,se.fit = TRUE, type="link")$fit
newdata$selink <- predict(avg_model, newdata, re.form=NA,se.fit = TRUE, type="link")$se.fit
newdata$realval <- (newdata$fitlink)
newdata$upper.ci <- (newdata$fitlink + 1.96*newdata$selink)
newdata$lower.ci <- (newdata$fitlink - 1.96*newdata$selink)
newdata$origlength<-newdata$s.DominantsAbundance * attr(E.bur$s.DominantsAbundance, 'scaled:scale') + attr(E.bur$s.DominantsAbundance, 'scaled:center')
newdata 

#graph
shannon <-ggplot(data= newdata, aes(x=origlength, y=realval), size = 5) + 
  labs(x = "", y = "SDI")+
  scale_x_continuous(expand = expansion(mult = c(0, 0))) + scale_y_continuous(expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#88419d") +
  geom_line() +
   theme_classic()+
  theme(axis.text.x=element_blank())+ # label left in blank to stack all graphs later
  theme(text = element_text(size = 18)) 


#Putting together the graphs (Figure 1)

EBpredictions <- plot_grid(shannon,spp, pielou, ncol=1, nrow = 3, labels = c("auto"), label_size = 14, align = "v")
EBpredictions
ggsave("EBpredictionsSept2025.png",EBpredictions,width=8,height = 10,dpi=300)


####============= L.praedator species=================#####

### 1) Evenness (Pielou)

#### a- Antswarm traits model

##### Beta

PielouLabA=glmmTMB(pielou.ev ~ s.velocity + s.biomass * s.total_area + (1|site),
                 data = L.pra, family=beta_family)
summary(PielouLabA)

#testing residuals
qqnorm(resid(PielouLabA),pch=16); qqline(resid(PielouLabA))

testDispersion(PielouLabA)
simulationOutput <- simulateResiduals(fittedModel = PielouLabA, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

#model selection
options(na.action= "na.fail")
mslist <- dredge(PielouLabA)
mslist
topmod<-get.models(mslist, subset=cumsum(weight) <= .95)
topmod # here you can see all of the models that make up 95% of the model weight
avgm <- model.avg(topmod)

summary(avgm)
confint(avgm, full = T) 

####b- Bird community model

##### Beta

L.pra$s.DominantsAbundance2 <- L.pra$s.DominantsAbundance^2
PielouLabB=glmmTMB(pielou.ev ~ s.DominantsAbundance*DominantBirds + s.DominantsAbundance2
                         +(1|site), data = L.pra, family=beta_family)
summary(PielouLabB)
confint(PielouLabB)


#testing residuals
qqnorm(resid(PielouLabB),pch=16); qqline(resid(PielouLabB))

testDispersion(PielouLabB)
simulationOutput <- simulateResiduals(fittedModel = PielouLabB, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

#model selection

options(na.action= "na.fail")
m1=glmmTMB(pielou.ev ~ s.DominantsAbundance*DominantBirds + s.DominantsAbundance2
                         +(1|site), data = L.pra, family=beta_family)
m2=glmmTMB(pielou.ev ~ s.DominantsAbundance+DominantBirds + s.DominantsAbundance2
                         +(1|site), data = L.pra, family=beta_family)
m3=glmmTMB(pielou.ev ~ s.DominantsAbundance + s.DominantsAbundance2
                         +(1|site), data = L.pra, family=beta_family)
m4=glmmTMB(pielou.ev ~ s.DominantsAbundance +(1|site), data = L.pra, family=beta_family)
m5=glmmTMB(pielou.ev ~ 1 +(1|site), data = L.pra, family=beta_family)
m6=glmmTMB(pielou.ev ~ s.DominantsAbundance*DominantBirds +(1|site), data = L.pra, family=beta_family)
m7=glmmTMB(pielou.ev ~ DominantBirds +(1|site), data = L.pra, family=beta_family)
m8=glmmTMB(pielou.ev ~ s.DominantsAbundance+DominantBirds +(1|site), data = L.pra, family=beta_family)

mslist<-list(m1,m2,m3,m4,m5,m6,m7,m8)

# Assign names 
names(mslist) <- paste0("m", 1:8)

# Get AICc table manually
aic_table <- model.sel(mslist)  # this returns a 'model.selection' object

# View AICc and weights
print(aic_table)

# Get cumulative sum of weights
cum_weights <- cumsum(aic_table$weight)

# Select models with cumulative weight ≤ 0.95
selected_models <- aic_table[cum_weights <= 0.95, ]

if (max(cum_weights[cum_weights <= 0.95]) < 0.95) {
  selected_models <- aic_table[1:(nrow(selected_models) + 1), ]
}

# Extract model names
selected_model_names <- rownames(selected_models)

# Subset model list
top_models <- mslist[selected_model_names]

# Perform model averaging on the top models
avg_model <- model.avg(top_models)

# Check the averaged model summary
summary(avg_model)

confint(avg_model, full = T)

### 2) Number of species

####a- Antswarm traits model

##### Poisson

PoissonSppLabA=glmmTMB(n_species ~ s.velocity + s.biomass * s.total_area +  offset(log(Total_SpeciesSite)) + (1|site),
                 data = L.pra, family=poisson)
summary(PoissonSppLabA)
confint(PoissonSppLabA)

#checking model residuals
qqnorm(resid(PoissonSppLabA),pch=16); qqline(resid(PoissonSppLabA))

testDispersion(PoissonSppLabA)
simulationOutput <- simulateResiduals(fittedModel = PoissonSppLabA, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

E1 <- resid(PoissonSppLabA, type = "pearson")
N  <- nrow(E.bur)
p  <- length(fixef(PoissonSppLabA)) + 1
Overdispersion <- sum(E1^2) / (N - p)
Overdispersion		


#model selection
options(na.action= "na.fail")
m1=glmmTMB(n_species ~ s.velocity + s.biomass * s.total_area +  offset(log(Total_SpeciesSite)) + (1|site),
                 data = L.pra, family=poisson)
m2=glmmTMB(n_species ~ s.velocity + s.biomass + s.total_area +  offset(log(Total_SpeciesSite)) + (1|site),
                 data = L.pra, family=poisson)
m3=glmmTMB(n_species ~ s.velocity + s.biomass  +  offset(log(Total_SpeciesSite)) + (1|site),
                 data = L.pra, family=poisson)
m4=glmmTMB(n_species ~ s.velocity + offset(log(Total_SpeciesSite)) + (1|site),
                 data = L.pra, family=poisson)
m5=glmmTMB(n_species ~ 1 + offset(log(Total_SpeciesSite)) + (1|site),
                 data = L.pra, family=poisson)
m6=glmmTMB(n_species ~ s.biomass * s.total_area +  offset(log(Total_SpeciesSite)) + (1|site),
                 data = L.pra, family=poisson)
m7=glmmTMB(n_species ~  s.biomass + s.total_area +  offset(log(Total_SpeciesSite)) + (1|site),
                 data = L.pra, family=poisson)
m8=glmmTMB(n_species ~ s.velocity + s.total_area +  offset(log(Total_SpeciesSite)) + (1|site),
                 data = L.pra, family=poisson)
m9=glmmTMB(n_species ~ s.total_area +  offset(log(Total_SpeciesSite)) + (1|site),
                 data = L.pra, family=poisson)
m10=glmmTMB(n_species ~  s.biomass  +  offset(log(Total_SpeciesSite)) + (1|site),data = L.pra, family=poisson)
  
mslist<-list(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10)

# Assign names 
names(mslist) <- paste0("m", 1:10)

# Get AICc table manually
aic_table <- model.sel(mslist)  

# View AICc and weights
print(aic_table)

# Get cumulative sum of weights
cum_weights <- cumsum(aic_table$weight)

# Select models with cumulative weight ≤ 0.95
selected_models <- aic_table[cum_weights <= 0.95, ]

if (max(cum_weights[cum_weights <= 0.95]) < 0.95) {
  selected_models <- aic_table[1:(nrow(selected_models) + 1), ]
}

# Extract model names
selected_model_names <- rownames(selected_models)

# Subset model list
top_models <- mslist[selected_model_names]

# Perform model averaging on the top models
avg_model <- model.avg(top_models)

# Check the averaged model summary
summary(avg_model)
confint(avg_model, full = T)

#### b- Bird community model

##### Poisson

PoissonSppLabB=glmmTMB(n_species ~ s.DominantsAbundance*DominantBirds + I(s.DominantsAbundance^2) +  offset(log(Total_SpeciesSite)) + (1|site),
                 data = L.pra, family=poisson)
summary(PoissonSppLabB)
confint(PoissonSppLabB)


#checking model residuals
qqnorm(resid(PoissonSppLabB),pch=16); qqline(resid(PoissonSppLabB))
testDispersion(PoissonSppLabB)
simulationOutput <- simulateResiduals(fittedModel = PoissonSppLabB, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

E1 <- resid(PoissonSppLabB, type = "pearson")
N  <- nrow(E.bur)
p  <- length(fixef(PoissonSppLabB)) + 1
Overdispersion <- sum(E1^2) / (N - p)
Overdispersion		

#model selection
options(na.action= "na.fail")
m1=glmmTMB(n_species ~ s.DominantsAbundance*DominantBirds + I(s.DominantsAbundance^2) +  offset(log(Total_SpeciesSite)) + (1|site),
                 data = L.pra, family=poisson)
m2=glmmTMB(n_species  ~ s.DominantsAbundance+DominantBirds + I(s.DominantsAbundance^2)+  offset(log(Total_SpeciesSite))
                         +(1|site), data = L.pra, family=poisson)
m3=glmmTMB(n_species ~ s.DominantsAbundance + I(s.DominantsAbundance^2)+  offset(log(Total_SpeciesSite))
                         +(1|site), data = L.pra, family=poisson)
m4=glmmTMB(n_species ~ s.DominantsAbundance+  offset(log(Total_SpeciesSite)) +(1|site), data = L.pra, family=poisson)
m5=glmmTMB(n_species ~ 1 +  offset(log(Total_SpeciesSite))+(1|site), data = L.pra, family=poisson)
m6=glmmTMB(n_species ~ s.DominantsAbundance*DominantBirds+ offset(log(Total_SpeciesSite)) +(1|site), data = L.pra, family=poisson)
m7=glmmTMB(n_species ~ DominantBirds +offset(log(Total_SpeciesSite))+(1|site), data = L.pra, family=poisson)
m8=glmmTMB(n_species~ s.DominantsAbundance+DominantBirds +  offset(log(Total_SpeciesSite))+(1|site), data = L.pra, family=poisson)

mslist<-list(m1,m2,m3,m4,m5,m6,m7,m8)

# Assign names 
names(mslist) <- paste0("m", 1:8)

# Get AICc table manually
aic_table <- model.sel(mslist)  # this returns a 'model.selection' object

# View AICc and weights
print(aic_table)

# Get cumulative sum of weights
cum_weights <- cumsum(aic_table$weight)

# Select models with cumulative weight ≤ 0.95
selected_models <- aic_table[cum_weights <= 0.95, ]

if (max(cum_weights[cum_weights <= 0.95]) < 0.95) {
  selected_models <- aic_table[1:(nrow(selected_models) + 1), ]
}

# Extract model names
selected_model_names <- rownames(selected_models)

# Subset model list
top_models <- mslist[selected_model_names]

# Perform model averaging on the top models
avg_model <- model.avg(top_models)

# Check the averaged model summary
summary(avg_model)
confint(avg_model, full = T)


### 3) Shannon

#### a- Antswarm traits model

##### Gaussian

GaussianShannonLabA=glmmTMB(shannon.di ~ s.velocity + s.biomass * s.total_area +  offset(log(max_ShannonSite)) + (1|site),
                 data = L.pra, family=gaussian)
summary(GaussianShannonLabA)
confint(GaussianShannonLabA)

#checking model residuals
qqnorm(resid(GaussianShannonLabA),pch=16); qqline(resid(GaussianShannonLabA))
testDispersion(GaussianShannonLabA)
simulationOutput <- simulateResiduals(fittedModel = GaussianShannonLabA, plot = F)
residuals(simulationOutput)
plot(simulationOutput)


#model selection
options(na.action= "na.fail")
m1=glmmTMB(shannon.di ~ s.velocity + s.biomass * s.total_area +  offset(max_ShannonSite) + (1|site),
                 data = L.pra, family=gaussian)
m2=glmmTMB(shannon.di ~ s.velocity + s.biomass + s.total_area +  offset(max_ShannonSite) + (1|site),
                 data = L.pra, family=gaussian)
m3=glmmTMB(shannon.di ~ s.velocity + s.biomass  +  offset(max_ShannonSite) + (1|site),
                 data = L.pra, family=gaussian)
m4=glmmTMB(shannon.di ~ s.velocity + offset(max_ShannonSite) + (1|site),
                 data = L.pra, family=gaussian)
m5=glmmTMB(shannon.di ~ 1 + offset(max_ShannonSite) + (1|site),
                 data = L.pra, family=gaussian)
m6=glmmTMB(shannon.di ~ s.biomass * s.total_area +  offset(max_ShannonSite) + (1|site),
                 data = L.pra, family=gaussian)
m7=glmmTMB(shannon.di ~  s.biomass + s.total_area +  offset(max_ShannonSite) + (1|site),
                 data = L.pra, family=gaussian)
m8=glmmTMB(shannon.di ~ s.velocity + s.total_area +  offset(max_ShannonSite) + (1|site),
                 data = L.pra, family=gaussian)
m9=glmmTMB(shannon.di ~ s.total_area +  offset(max_ShannonSite) + (1|site),
                 data = L.pra, family=gaussian)
m10=glmmTMB(shannon.di ~  s.biomass  +  offset(max_ShannonSite) + (1|site),data = L.pra, family=gaussian)
  
mslist<-list(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10)

# Assign names 
names(mslist) <- paste0("m", 1:10)

# Get AICc table manually
aic_table <- model.sel(mslist)  # this returns a 'model.selection' object

# View AICc and weights
print(aic_table)

# Get cumulative sum of weights
cum_weights <- cumsum(aic_table$weight)

# Select models with cumulative weight ≤ 0.95
selected_models <- aic_table[cum_weights <= 0.95, ]

if (max(cum_weights[cum_weights <= 0.95]) < 0.95) {
  selected_models <- aic_table[1:(nrow(selected_models) + 1), ]
}

# Extract model names
selected_model_names <- rownames(selected_models)

# Subset model list
top_models <- mslist[selected_model_names]

# Perform model averaging on the top models
avg_model <- model.avg(top_models)

# Check the averaged model summary
summary(avg_model)
confint(avg_model, full = T)

#### b- Bird community model

##### Gaussian
GaussianShannonLabB=glmmTMB(shannon.di ~ s.DominantsAbundance*DominantBirds + I(s.DominantsAbundance^2) +  offset(log(max_ShannonSite)) + (1|site),
                 data = L.pra, family=gaussian)
summary(GaussianShannonLabB)
confint(GaussianShannonLabB)

#checking model residuals
qqnorm(resid(GaussianShannonLabB),pch=16); qqline(resid(GaussianShannonLabB))

testDispersion(GaussianShannonLabB)
simulationOutput <- simulateResiduals(fittedModel = GaussianShannonLabB, plot = F)
residuals(simulationOutput)
plot(simulationOutput)


#model selection
options(na.action= "na.fail")

m1=glmmTMB(shannon.di ~ s.DominantsAbundance*DominantBirds + I(s.DominantsAbundance^2) +  offset(max_ShannonSite) + (1|site),
                 data = L.pra, family=gaussian)
m2=glmmTMB(shannon.di  ~ s.DominantsAbundance+DominantBirds + I(s.DominantsAbundance^2)+  offset(max_ShannonSite)
                         +(1|site), data = L.pra, family=gaussian)
m3=glmmTMB(shannon.di ~ s.DominantsAbundance + I(s.DominantsAbundance^2)+  offset(max_ShannonSite)
                         +(1|site), data = L.pra, family=gaussian)
m4=glmmTMB(shannon.di ~ s.DominantsAbundance+  offset(max_ShannonSite) +(1|site), data = L.pra, family=gaussian)
m5=glmmTMB(shannon.di ~ 1 +  offset(max_ShannonSite)+(1|site), data = L.pra, family=gaussian)
m6=glmmTMB(shannon.di ~ s.DominantsAbundance*DominantBirds+ offset(max_ShannonSite) +(1|site), data = L.pra, family=gaussian)
m7=glmmTMB(shannon.di ~ DominantBirds +offset(max_ShannonSite)+(1|site), data = L.pra, family=gaussian)
m8=glmmTMB(shannon.di~ s.DominantsAbundance+DominantBirds +  offset(max_ShannonSite)+(1|site), data = L.pra, family=gaussian)

mslist<-list(m1,m2,m3,m4,m5,m6,m7,m8)

# Assign names 
names(mslist) <- paste0("m", 1:8)

# Get AICc table manually
aic_table <- model.sel(mslist)  # this returns a 'model.selection' object

# View AICc and weights
print(aic_table)

# Get cumulative sum of weights
cum_weights <- cumsum(aic_table$weight)

# Select models with cumulative weight ≤ 0.95
selected_models <- aic_table[cum_weights <= 0.95, ]

if (max(cum_weights[cum_weights <= 0.95]) < 0.95) {
  selected_models <- aic_table[1:(nrow(selected_models) + 1), ]
}

# Extract model names
selected_model_names <- rownames(selected_models)

# Subset model list
top_models <- mslist[selected_model_names]

# Perform model averaging on the top models
avg_model <- model.avg(top_models)

# Check the averaged model summary
summary(avg_model)
confint(avg_model, full = T)


##====== Do antswarm traits influence the number of dominants?=======##

### E. burchellii

#how do swarm traits influence number of dominants
#random effects selection
m1=glmmTMB(DominantsAbundance ~ s.velocity + s.biomass * s.total_area + (1|swarm_id),
                 data = E.bur, family=gaussian)
m2=glmmTMB(DominantsAbundance ~ s.velocity + s.biomass * s.total_area + (1|site),
                 data = E.bur, family=gaussian)
m3=glmmTMB(DominantsAbundance ~ s.velocity + s.biomass * s.total_area + (1|site) + (1|swarm_id),
                 data = E.bur, family=gaussian)
m4=glmmTMB(DominantsAbundance ~ s.velocity + s.biomass * s.total_area + (1|site/swarm_id),
                 data = E.bur, family=gaussian)

summary(m1)
summary(m2)
summary(m3)
summary(m4)

#In this case both random effects, site and swarm id should be considered in the model

# global model
gmDominanceEciton=glmmTMB(DominantsAbundance ~ s.velocity + s.biomass * s.total_area + (1|site/swarm_id),
                 data = E.bur, family=gaussian)
summary(gmDominanceEciton)

#testing residuals
qqnorm(resid(gmDominanceEciton),pch=16); qqline(resid(gmDominanceEciton))
testDispersion(gmDominanceEciton)
simulationOutput <- simulateResiduals(fittedModel = gmDominanceEciton, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

#model selection
options(na.action= "na.fail")
mslist <- dredge(gmDominanceEciton)
mslist
topmod<-get.models(mslist, subset=cumsum(weight) <= .95)
topmod # here you can see all of the models that make up 95% of the model weight
avgm <- model.avg(topmod)

summary(avgm)
confint(avgm, full = T) 

# Create the scaled variable correctly
s.total_area <- seq(
  from = min(E.bur$s.total_area),
  to = quantile(E.bur$s.total_area, 0.95),
  by = 0.1
)

#generate predictions
s.total_area=seq(from=min(E.bur$s.total_area),to=quantile(E.bur$s.total_area,c(0.95)),by=0.1)
s.biomass= mean(E.bur$s.biomass)
s.velocity= mean(E.bur$s.velocity)

newdata <- data.frame(expand.grid(site=unique(E.bur$site),swarm_id=unique(E.bur$swarm_id), s.biomass=s.biomass, s.velocity=s.velocity, s.total_area=s.total_area))

newdata$fitlink <- predict(avgm, newdata,re.form=NA,se.fit = TRUE, type="link",allow.new.levels=TRUE)$fit
newdata$selink <- predict(avgm, newdata, re.form=NA,se.fit = TRUE, type="link",allow.new.levels=TRUE)$se.fit
newdata$realval <- (newdata$fitlink)
newdata$upper.ci <- (newdata$fitlink + 1.96*newdata$selink)
newdata$lower.ci <- (newdata$fitlink - 1.96*newdata$selink)
newdata$origlength<-newdata$s.total_area * attr(E.bur$s.total_area, 'scaled:scale') + attr(E.bur$s.total_area, 'scaled:center')
newdata 

#graph
dom_plotting_data <- subset(E.bur, E.bur$total_area < quantile(E.bur$total_area, 0.95))

domabund <-ggplot(data= newdata, aes(x=origlength, y=realval), size = 5) + 
  labs(x = expression("Swarm area (m"^2*")"), y = "Number of dominant individuals")+
  scale_x_continuous(expand = expansion(mult = c(0, 0))) + scale_y_continuous(expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#88419d") +
  geom_line() +
  theme_classic()+ 
  theme(text = element_text(size = 18))
domabund
ggsave("domabundSept2025.png", plot = domabund, width = 8, height = 6, dpi = 150)


### L. praedator model

hist(L.pra$DominantsAbundance)

# We have so many zeros we are going to do a binomial for presence/absence of a dominant birds in the L. praedator swarms.

L.pra1 <- L.pra %>% 
  mutate(DominantsAbundance= replace(DominantsAbundance, DominantsAbundance>0, 1))

gmDominanceLabidus=glmmTMB(DominantsAbundance ~ s.velocity + s.biomass * s.total_area + (1|site),
                 data = L.pra1, family=binomial)
summary(gmDominanceLabidus)

#testing residuals
qqnorm(resid(gmDominanceLabidus),pch=16); qqline(resid(gmDominanceLabidus))

testDispersion(gmDominanceLabidus)
simulationOutput <- simulateResiduals(fittedModel = gmDominanceLabidus, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

#model selection
options(na.action= "na.fail")
mslist <- dredge(gmDominanceLabidus)
mslist
topmod<-get.models(mslist, subset=cumsum(weight) <= .95)
topmod # here you can see all of the models that make up 95% of the model weight
avgm <- model.avg(topmod)

summary(avgm)
confint(avgm, full = T) 



###=========== What influences the number of non-domiants?================## 

#E. burchellii
#First, have to add one column with the number of Non-dominant birds. 

EBdom <- E.bur %>% 
 mutate(NonDom= round(n_individuos-DominantsAbundance)) 

#test for different random effects

m1 <- glmer.nb(NonDom ~ s.DominantsAbundance + I(s.DominantsAbundance^2)+ s.velocity + s.biomass * s.total_area + (1|site),
                 data = EBdom)
m2 <- glmer.nb(NonDom ~ s.DominantsAbundance + I(s.DominantsAbundance^2)+ s.velocity + s.biomass * s.total_area + (1|swarm_id),
                 data = EBdom)
m3 <- glmer.nb(NonDom ~ s.DominantsAbundance + I(s.DominantsAbundance^2)+ s.velocity + s.biomass * s.total_area + (1|site) + (1|swarm_id),
                 data = EBdom)
m4 <- glmer.nb(NonDom ~ s.DominantsAbundance + I(s.DominantsAbundance^2)+ s.velocity + s.biomass * s.total_area + (1|site/swarm_id),
                 data = EBdom)

summary(m1)
summary(m2)
summary(m3)
summary(m4)
#convergence issues except for one with site

#final global model
NonDom<- glmmTMB(NonDom ~ s.DominantsAbundance + I(s.DominantsAbundance^2)+ s.velocity + s.biomass * s.total_area + (1|site),
                 family = nbinom2, data=EBdom)
summary(NonDom)

#test residuals
qqnorm(resid(NonDom),pch=16); qqline(resid(NonDom))
testDispersion(NonDom)
simulationOutput <- simulateResiduals(fittedModel = NonDom, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

#model selection
options(na.action= "na.fail")
mslist <- dredge(NonDom, subset = dc(s.DominantsAbundance, I(s.DominantsAbundance^2)))
mslist
topmod<-get.models(mslist, subset=cumsum(weight) <= .95)
topmod # here you can see all of the models that make up 95% of the model weight
avgm <- model.avg(topmod)

summary(avgm)
confint(avgm, full = T) 

#predictions
quantile(EBdom$s.DominantsAbundance,c(0.1,0.5,0.90,0.95)) 

scale_center <- attr(EBdom$s.DominantsAbundance, "scaled:center")
scale_scale <- attr(EBdom$s.DominantsAbundance, "scaled:scale")

s.DominantsAbundance=seq(from=min(EBdom$s.DominantsAbundance),to=quantile(EBdom$s.DominantsAbundance,c(0.95)) ,by=0.1)
s.biomass= mean(EBdom$s.biomass)
s.velocity= mean(EBdom$s.velocity)
s.total_area= mean(EBdom$s.total_area)

newdata <- data.frame(expand.grid(site=unique(EBdom$site), s.DominantsAbundance=s.DominantsAbundance,s.biomass=s.biomass, s.velocity=s.velocity, s.total_area=s.total_area))

newdata$fitlink <- predict(avgm, newdata,re.form=NA,se.fit = TRUE, type="link")$fit
newdata$selink <- predict(avgm, newdata, re.form=NA,se.fit = TRUE, type="link")$se.fit
newdata$realval <- exp(newdata$fitlink)
newdata$upper.ci <- exp(newdata$fitlink + 1.96*newdata$selink)
newdata$lower.ci <- exp(newdata$fitlink - 1.96*newdata$selink)
newdata$origlength<-newdata$s.DominantsAbundance * attr(E.bur$s.DominantsAbundance, 'scaled:scale') + attr(E.bur$s.DominantsAbundance, 'scaled:center')
newdata 

#graph

nondomabund <-ggplot(data= newdata, aes(x=origlength, y=realval), size = 5) + 
  labs(x = expression("Number of dominant individuals"), y = "Number of non-dominant individuals")+
  scale_x_continuous(expand = expansion(mult = c(0, 0))) + scale_y_continuous(expand = expansion(mult = c(0, 0)))+
  geom_ribbon(aes(ymin = lower.ci, ymax = upper.ci), alpha = .45, fill="#88419d") +
  geom_line() +
  theme_classic()+ 
  theme(text = element_text(size = 18))
nondomabund

dompredictions <- plot_grid(nondomabund,domabund, ncol=2, nrow = 1, labels = c("auto"), label_size = 12, align = "hv")
dompredictions
ggsave("dompredictionsSept2025test.png",dompredictions,width=9,height = 6,dpi=300)
#above is Fig 2 in the paper


## L. praedator

LPdom <- L.pra %>% 
  mutate(NonDom= round(n_individuos-DominantsAbundance)) 

# antswarm traits and non-dominant abundance

#find random effect structure

m1 <- glmmTMB(NonDom ~ s.velocity + s.biomass * s.total_area + (1|site),family = nbinom2,
                 data = LPdom)
m2 <- glmmTMB(NonDom ~ s.velocity + s.biomass * s.total_area + (1|swarm_id),family = nbinom2,
                 data = LPdom)
m3 <- glmmTMB(NonDom ~ s.velocity + s.biomass * s.total_area + (1|swarm_id)+ (1|site),family = nbinom2,
                 data = LPdom)
m4 <- glmmTMB(NonDom ~ s.velocity + s.biomass * s.total_area + (1|site/swarm_id),family = nbinom2,
                 data = LPdom)
summary(m1)
summary(m2)
summary(m3)
summary(m4)


#final global model
NonDomlab<- glmmTMB(NonDom ~ s.velocity + s.biomass * s.total_area + (1|site),family = nbinom2, data=LPdom)
summary(NonDomlab)

#residuals
qqnorm(resid(NonDomlab),pch=16); qqline(resid(NonDomlab))
testDispersion(NonDomlab)
simulationOutput <- simulateResiduals(fittedModel = NonDomlab, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

#model selection
options(na.action= "na.fail")
mslist <- dredge(NonDomlab)
mslist
topmod<-get.models(mslist, subset=cumsum(weight) <= .95)
topmod # here you can see all of the models that make up 95% of the model weight
avgm <- model.avg(topmod)

summary(avgm)
confint(avgm, full = T) 


# Now with the dominant abundance effects on number of nondominants

#test different random effects
m1 <- glmmTMB(NonDom ~ s.DominantsAbundance + I(s.DominantsAbundance^2) + (1|site),family = nbinom2,
                 data = LPdom)
m2 <- glmmTMB(NonDom ~ s.DominantsAbundance + I(s.DominantsAbundance^2) + (1|swarm_id),family = nbinom2,
                 data = LPdom)
m3 <- glmmTMB(NonDom ~ s.DominantsAbundance + I(s.DominantsAbundance^2)+ (1|swarm_id)+ (1|site),family = nbinom2, data = LPdom)

m4 <- glmmTMB(NonDom ~ s.DominantsAbundance + I(s.DominantsAbundance^2) + (1|site/swarm_id),family = nbinom2,
                 data = LPdom)

summary(m1)
summary(m2)
summary(m3)
summary(m4)

#final global model
NonDomlabB<- glmmTMB(NonDom ~ s.DominantsAbundance + I(s.DominantsAbundance^2) + (1|site),family = nbinom2, data=LPdom)
summary(NonDomlabB)

#model residuals
qqnorm(resid(NonDomlabB),pch=16); qqline(resid(NonDomlabB))

testDispersion(NonDomlabB)
simulationOutput <- simulateResiduals(fittedModel = NonDomlabB, plot = F)
residuals(simulationOutput)
plot(simulationOutput)

#model selection
options(na.action= "na.fail")
m1<- glmmTMB(NonDom ~ s.DominantsAbundance + I(s.DominantsAbundance^2) + (1|site),family = nbinom2, data=LPdom)
m2<- glmmTMB(NonDom ~ s.DominantsAbundance  + (1|site),family = nbinom2, data=LPdom)
m3<- glmmTMB(NonDom ~ 1  + (1|site),family = nbinom2, data=LPdom)

mslist<-list(m1,m2,m3)

# Assign names 
names(mslist) <- paste0("m", 1:3)

# Get AICc table manually
aic_table <- model.sel(mslist)  # this returns a 'model.selection' object

# View AICc and weights
print(aic_table)

# Get cumulative sum of weights
cum_weights <- cumsum(aic_table$weight)

# Select models with cumulative weight ≤ 0.95
selected_models <- aic_table[cum_weights <= 0.95, ]

if (max(cum_weights[cum_weights <= 0.95]) < 0.95) {
  selected_models <- aic_table[1:(nrow(selected_models) + 1), ]
}

# Extract model names
selected_model_names <- rownames(selected_models)

# Subset model list
top_models <- mslist[selected_model_names]

# Perform model averaging on the top models
avg_model <- model.avg(top_models)

# Check the averaged model summary
summary(avg_model)
confint(avg_model, full = T)


























